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Abstract. The time-dependent surface flux (t-SURFF) method is introduced for 
computing of strong-field infrared photo-ionization spectra of atoms by numerically 
solving the time-dependent Schrodinger equation on minimal simulation volumes. The 
volumes only need to accommodate the electron quiver motion and the relevant range 
of the atomic binding potential. Spectra are computed from the electron flux through a 
surface, beyond which the outgoing flux is absorbed by infinite range exterior complex 
scaling (irECS). Highly accurate infrared photo-electron spectra are calculated in single 
active electron approximation and compared to literature results. Detailed numerical 
evidence for performance and accuracy is given. Extensions to multi-electron systems 
and double ionization are discussed. 
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1. Introduction 

In a broad range of recent experiments, strong infrared laser pulses, often combined with 
high harmonic pulses, are used to study the electronic dynamics of atoms and molecules 
on the natural time scale of valence electron motion of < 1 fs. Basic mechanisms 
of the IR-electron interaction are well understood within the simple semi-classical re- 
collision model [1], but for a more detailed understanding numerical simulations must 
be employed. This is due to the fundamentally non-perturbative interaction of near IR 
fields with valence electrons at intensities of > 10 14 W/cm 2 . Even for the simplest 
single-electron models the simulation remains challenging, especially when accurate 
photo-electron momentum spectra are required, as, e.g., for re-collision imaging [2, 3, 4]. 
When two-electron processes are involved, one quickly reaches the limits of present day 
computer resources [5, 6]. 

The surprising difficulty in simulating a seemingly simple process like ionization by 
a dipole field is due to the presence of vastly different length- and time-scales: first, even 
though the laser pulses can be as short as a single optical cycle, at the Ti:Sapphire wave 
length of A = 800 nm this still corresponds to a FWHM duration of > 2.5 fs or about 
110 atomic units (a.u., h = e 2 = m e = 1). During that time, a photo-electron with 
an energy of 13 eV ~ 1/2 a.u. moves to a distance of ~ 110 Bohr, which sets a lower 
limit for the required box-size, if reflections from box boundaries are to be avoided. 
In practice, higher energies and longer pulse durations including the rise and fall of 
the pulses are of interest, leading to simulation volumes with diameters of thousands 
of atomic units. At the same time, photo-electron spectra are broad, extending to at 
least 2 U p for the "direct" photo-electrons and further up to about 10 U p for the "re- 
scattered" electrons. Re-scattered electrons are those that after ionization re-directed to 
the nucleus by the laser field, where they absorb more photons in an inelastic scattering 
process. The ponderomotive potential U p = I/(Au 2 ) grows linearly with laser intensity 
/ and quadratically with wavelength. At the moderate intensity / = 10 lA W/cm 2 and 
A = 800 nm it is U p = 0.22 a.u. ~ 6eV. The re-scattering momentum energy cutoff 
of 10U P corresponds to a photo-electron momentum of 2.2 a.u. For representation of 
such momenta on a spatial grid, we need grid spacings of at least Ax < 2n/2 a.u., for 
accurate results usually significantly more than this. This leaves us with thousands of 
grid points in each spatial direction even for moderate laser parameters. The situation 
quickly worsens at higher intensities and longer wavelengths. 

This general requirement on discretization cannot be overcome by any specific 
representation of the wave function: speaking in terms of classical mechanics, we must 
represent the phase space that is covered by the electrons, which involves a certain 
range of momenta and positions. If we have no additional knowledge of the structure of 
solution, the number of discretization points we need is the phase space volume divided 
by the Planck constant h. In some cases like, for example, single-photon ionization, 
we can exploit the fact that at long distances the solution covers only a very narrow 
range of momenta and only the spatially well-localized initial bound state requires a 
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broader range of momenta: the phase-space volume remains small, simple models like 
perturbation theory allow reproducing the physics. We have no such simplifying physical 
insight for strong-field IR photo-ionization. 

The lower limit for the number of discretization points for the complete wave 
function can be approached by different strategies: the choice of velocity gauge [7], 
working in the Kramers- Henneberger frame [8] or in momentum space [9], by variable 
grid spacings, or by expanding into time- dependent basis functions [10]. A promising 
strategy is to follow the solution in time [1 1] . 

Alternatively, we can abandon the attempt of representing the complete wave- 
function and instead use absorbing boundaries and extract momenta at finite distances. 
The time-dependent surface flux method (t-SURFF) introduced here is such an 
approach. After its mathematical derivation, numerical implementation is briefly 
discussed. Angle-resolved photo-electron spectra are presented using between 75 and 
200 radial discretization points for truncated and full Coulomb potentials. We discuss 
accuracies and demonstrate the efficiency of t-SURFF by comparison with recent 
literature. Finally, possible extensions to few-electron systems and double photo- 
electron spectra are outlined. 



2. The t-SURFF method 



Scattering measurements and theory are both based on the plausible idea that 
interactions are limited to finite ranges in space and time and that at large times > T 
and large distances > R the time-evolution of the scattering particle is that of free 
motion: 

V(f,t) ~ J dk (3) exp(-itk 2 /2)b{k)x%{f) fort >T,\f\> R, (1) 
where Xkif) — (27r)~ 3//2 exp(ik ■ r) are ^-normalized plane waves. The measured 

— * 

momentum spectrum is proportional to the square of the spectral amplitudes b(k) 

a(k) oc \b\k)\. (2) 
For Hamiltonians that are time-independent beyond a certain time T, we can readily 

— * 

obtain the spectral amplitudes b(k) by decomposing the wave function \P(f, T) into its 
spectral components 

b(k) = (^(T))eM-^k 2 /2), (3) 

where the scattering solutions 



2 



have the asymptotic behavior 

^(r) ~ Xfe(r) for M -> 00 ■ ( 5 ) 

For computing b(k), one needs to (i) propagate a solution ^(f,t) until time T, (ii) 
obtain the scattering solutions ip^. Unfortunately, both of these tasks are non-trivial 
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in all but the simplest cases. As discussed above, solving the TDSE with IR fields is 
numerically challenging because of the large box sizes needed. "Obtaining the scattering 
solution" amounts to outright solving a time-independent scattering problem for H(T), 
but analytic scattering wave functions are only known for simple model potentials and 
for the Coulomb potential. 

Rather than letting the system evolve and analyzing it at the end of the evolution, 
we can record the particle flux leaving a finite volume as the systems evolves. Such a 
procedure neither requires stationary scattering solutions nor the complete wave function 
in an asymptotically large volume. It is practical, if we either know the further time- 
evolution of the system outside the finite volume in analytic form or if it can be obtained 
with little numerical effort. This type of methods has been applied for reactive scattering 
with time-independent Hamiltonians [12, 13], where obtaining scattering solutions would 
be tantamount to solving the complete stationary scattering problem, a daunting task 
for few-body systems. 

Photo-ionization cannot be computed by existing surface flux methods, as the dipole 
interaction is non-local and the external field modifies the particle energies everywhere, 
in particular also after the particle has left the finite simulation volume. To handle 
this, we have developed the time-dependent surface flux (t-SURFF) method, that will 
be derived in the following. A preliminary version of the method was published in [14]. 

Let us choose a surface radius R c large enough that the particle motion can be 
considered a free motion and that all occupied bound states of the system have negligible 
particle density at \f\ > R c . Let us further pick a sufficiently large time T such that all 
particles that will ever reach our detector with energy k 2 /2 > are outside the finite 
volume \r\ < R c . At that time, the wave function has split into bound and asymptotic 
parts 

*(r,T) = V h (r,T) + y s (r,T) (6) 

with 

# b (r,T)«0 for|rl>i? c (7) 

* s (f,T) := J dk {3) exp(-iTk 2 /2)b(k)^(r) « for \f\ < R c . (8) 

The approximate sign in (8) refers to the fact that very low energy particles k 2 /2 ~ 
may not have left the finite volume at time T. It follows that, the lower the energy, the 
larger T will be required for the splitting to hold. Although each individual \ k extends 
over the complete space, the scattering wave-packet is localized outside R c up to a small 
error that quickly decays with growing T. The scattering amplitudes can be obtained 
as 

e* Tp l%{k) = (^|^s(T)) » (i> k \0(R c )\y s (T)) = ( Xk \e(R c )\* s (T)). (9) 
Here we introduced the notation 

(^ k \6(R c )\%(T)} := [ rf( 3 V^(f)M/ s (f,T) (10) 

J\P\>Rc 
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The substitution of the scattering solution with a plane wave — > \k m the last step 
uses the asymptotic behavior (5). It is exact, when all interactions vanish beyond R c . 

For converting the above matrix element to a time-integral over surface values, we 
must assume that we know the time-evolution of the particle after it has passed through 
the surface. For that we assume that there is a "channel Hamiltonian" H c (t) such that 

H c (t) = H{t) for \f] > R c andVt. (11) 

In case of a short range potential V(r) = for \r\ > R c this is exactly fulfilled by the 
Hamiltonian for the free motion in the laser field 

H e (t) = h-iv-A(t)] 2 , 



(12) 



where A(t) = — f^S^dt' for the dipole field £{t). The Volkov solutions 



x -(f) = (27r)- 3 /V^ fc 'V fc - r ; 



^{k,t)= l -j\r[k-A{t)f (13) 



solve the TDSE 

d 
%— 
dt 



Xi(t))=H c (t)\XiV)). 



(14) 



We can now write 

( Xk (T)\0(R c )\* s (T)) 



dtf t ( Xk (t)\9(R c )\* s (t)) 



dt( Xk (t)\H c (t)6(R c ) - 9(R c )H(t)\V s (t)) 



= i 



T 



dt( X k(t)\ 



--A + iA(t)-V,e(R c 



(15) 



The commutator vanishes everywhere except on the surface \r\ = R c . Assuming linear 
polarization in ^-direction A{t) = (0,0, A(t)), it can be written in polar coordinates 
(r, 9, 0) as 



--A + iA(t)d z ,6(R c ) 



1 



- -r- 2 d r r 2 5(r - R c ) - ^5(r - R c )d r - iA(t) cos 95(r - R c ) (16) 

With this, the volume integral over the space covered by the solution at time T has 
been converted to a time-integral up to T and a surface integral over \r\ = R c . 

We would like to remark that without time- dependence we can make one more step, 
as then the Volkov phase reduces to Q(k,t) = tk 2 /2 and time-integration turns into the 
time-energy Fourier-transform of the surface integral, which connects to the well-known 
results of, e.g., Ref. [13]. 
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3. Finite element discretization and irECS 

For the efficient use of t-SURFF, a reliable mechanism for truncating the solution outside 
the finite volume without generating reflections or other artefacts is needed. Commonly 
complex absorbing potentials are employed for that purpose (for a recent review, see 
[15]). We use infinite range exterior scaling (irECS) introduced and discussed in detail 
in reference [16]. Here we only briefly summarize the procedure. 

Exterior complex scaling (ECS) consists in making an analytic continuation of 
the Hamiltonian by rotating the coordinates into the upper complex plane beyond the 
"scaling radius" Rq\ 

f 6Ro = { f for \f\ < i? ^ [i? + e ie (\r\ - R )] for|f| > R (17) 

The resulting complex scaled Hamiltonian Hg Ro can be used in a complex scaled time- 
dependent Schrodinger equation 

t^eR Q (t) = H eRo (t)^ eRi) (t). (18) 

It was observed in [16] that for the velocity form of the TDSE, in the unsealed region 
\r\ < R exact and complex scaled solution agree: ^^(r, t) = ^ v (f,t). Mathematical 
proof is absent, but numerical agreement can be pushed to machine precision with 
relative errors ~ 10~ 14 . 

These high accuracies can be reached with little effort by infinite range ECS 
(irECS), where very few ~ 20 discretization coefficients are needed at radii > Rq. It 
consists in using an expansion into spherical harmonics and discretize the radial parts 
by high order finite elements. As the last element, irECS uses the infinite interval 
[i?o, oo) and functions of the form L n (2ar) exp(— ar) for its discretization, where L n are 
Laguerre polynomials. The idea is that low momentum / long de-Broglie wave length 
electrons need rather large absorption ranges to become absorbed, but that the solution 
at large distances is rather smooth. The more oscillatory short wave length content 
of the wave function becomes absorbed within a few oscillations. The exponentially 
damped functions turned out be very efficient in emulating that behavior: in most 
cases, about 20 functions at \r\ > Rq are sufficient for complete absorption. Typical 
finite element orders used are 15 to 25. The comparatively high order further enhances 
efficiency of irECS. 

Although in irECS there is no strict box boundary, the number of discretization 
points and the phase space volume that can be represented remain finite. There is a 
correlation between position discretization and momentum discretization in irECS: large 
electron momenta can only be represented in the unsealed region and at the beginning 
of the scaled region. At long distances into the scaled region, only gentle oscillations and 
therefore low momenta can be represented until the wave function dies out exponentially. 

As the finite element basis is local, the discretization coefficients are approximately 
associated with positions in space. In fact, we can establish a one-to-one relation between 
the M functions belonging an element and the function values at M different points 
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within that element. In that sense we will refer to the discretization coefficients as 
"discretization points" , emphasizing the locality of the finite element functions. 

The exact choice of the scaling parameters 9, Ro and a is not critical. In Ref. [16] 
we found little variation of accuracy for values in the intervals 9 ~ [0.5, 0.9] and 
a ~ [0.2, 0.6]. Variations of the results with R , as a rule, reflect the spatial discretization 
error of the calculation. We found this behavior confirmed also for calculation of the 
photo-electron spectra presented here. The plots shown below were all calculated with 
9 = a = 0.5. 

4. Photo-electron momentum spectra for a short range potential 

We solve the time-dependent Schrodinger equation in velocity gauge 

l Jt^ vm = ~ Ait)]2 + V{f) } (19) 

where we first use the short-range "Coulomb" potential 

V(r) = jc [-1/r -r 2 /(2R 3 ) + 3/(2/?)] forr < R Oforr > R. (20) 

With R = 20 and an effective charge c = 1.1664 the ground state energy is —0.5. The 
laser pulse is linearly polarized in z-direction with the vector potential 

A z (t) = -cos 2 A sin(tu). (21) 

We choose parameters u = 0.057 and So = 0.0755 corresponding to laser wave length 
of 800 nm and peak intensity 2 x 10 U W / cm 2 . T is the full width at half maximum 
(FWHM) of the vector potential, total pulse duration is 2T. 

Figure 1 shows the total and partial wave photo-electron spectra for potential range 
R = 20 and T = 5 optical cycles. At these parameters, more than 90% of the electrons 
get detached. For accuracy < 1% up to energies of 10 U p ~ 120 eV we need L max = 30 
partial waves with only 90 radial discretization points. We define the error relative to 
an accurate reference spectrum a re f as 

^ - J&i.wgL) - {am < E '-mC dE '° (E,) - (22) 

Including in the denominator the average over the interval [E — 5E, E + SE] suppresses 
spurious spikes in the error due to near-zeros of the spectrum. We choose 5E = 
0.05a.u. ~ 1.5 eV, which at 800 nm corresponds to averaging over about 2 photo-electron 
peaks. 

In the unsealed region we use 60 points, 30 points are located in \r\ > R - The 
accuracy estimate shown in Fig. 1 is obtained by comparing to a fully converged 
calculation. When we increase the number of points to 180, the error drops to < 10 -3 
The increase of relative errors with energy can be attributed to the decrease of the signal: 
note that from to 10U P the spectrum drops by more than 5 orders of magnitude. 

As a consistency check, we find that spectra up to 10 U p computed with largely 
different surface radii R c = 21 and 29 coincide within relative accuracies of better than 
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Figure 1. Upper panels: total and partial wave photo-electron spectra obtained 
for the smoothly truncated Coulomb potential Eq. (20). The lower panels show 
the relative errors for each spectrum according to Eq. (22) of a calculation with 
with only 90 discretization points per angular momentum comparing to a fully 
converged calculation. Pulse parameters: A = 800 nm, T = 5 optical cycles, intensity 
= 2 x 10 14 W/cm 2 . 



< 10~ 3 . As beyond R c we assume the exact Volkov propagation, this demonstrates that 
the solution of the TDSE and the integrals (15) are correct and that also reflections are 
suppressed on at least that level of accuracy. 

The calculated spectra are independent of the complex scaling parameters: over 
the ranges Rq G [20, 30] and 9 = [0.4, 0.7] results vary by less than 10~ 3 . The lower limit 
for R Q is not dictated by complex scaling: rather, as we want to obtain exact results, we 
must pick up the exact wave-function outside the range of the potential R c > R = 20, 
and therefore also Rq > R c > R. Already with these parameters the quiver amplitude, 
i.e. the excursion of free electrons in the laser field Sq/u> 2 ~ 23 reaches beyond Rq and 
into the complex scaled region. This confirms an earlier observation that the dynamics 
is correctly reproduced also in the complex scaled region [16]. 

For correct electron spectra, the effective box size of the combined unsealed and 
scaled regions must be large enough to accommodate the quiver motion. To study this 
further, we use a somewhat shorter effective range of R = 15 and choose Rq = R c = R. 
From Fig. 2 we see that with only 45 discretization point in the unsealed region r < R 
1% accuracy is reached at the same laser parameters as before. With 30 points for 
absorptions we have a total of 75 points. Note that the quiver radius of ~ 23 a.u. now 
reaches rather deep into scaled region but still fits into the total box. 

Keeping the intensity fixed, a longer wave length of 1000 nm leads to a larger quiver 
radius of ~ 36 a.u.. We expect to need a factor (1000/800) 2 ~ 1.6 more discretization 
points. Indeed, we can reach < 1% accuracy with total of 120 points in this case (Fig. 2). 
All additional coefficients and at least half of the quiver motion now are located in the 
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Photoelectron energy (eV) Photoelectron energy (eV) 

Figure 2. Required box sizes depend on the quiver amplitude. Left panels: spectrum 
(top) for a 5 cycle FWHM pulse at 800 nm and £ = 0.0755 (intensity 2 x 10 u W/cm 2 ). 
Accuracies (bottom) are calculated with 75 (black line) and 100 (red line) discretization 
points. Errors < 1% relative to an accurate reference calculation are reached with 75 
points. At 100 points the error is < 10~ 3 . Right panels: spectrum (top) and accuracies 
(bottom) for 5 cycles FWHM at 1000 nm and the same intensity. Because of the larger 
quiver amplitude, a ~ 60% larger box with 120 points is needed for < 1% accurate 
results. When using a smaller box with 120 points at a 30% higher density errors 
increase to near 10% (red line). 

scaled region. Note that also U p and with it the peak momentum grows with wave 
length: for describing the energies > 120 eV we would also need to increase the density 
of points by 25 %. However, just increasing the number of points without increasing 
the box size gives incorrect results: it appears that we need to accommodate the full 
quiver motion up to ~ 36 a. u. in the simulation box. 

5. Photo-electron momentum spectra for the hydrogen atom 

As always in scattering problems, the long-range nature of the Coulomb potential 
introduces extra mathematical and practical complications. Considering t-SURFF, 
there is no surface radius R c such that the Volkov solutions become exact. In addition, 
the Rydberg bound states extend to arbitrarily large distances. We will discuss below 
how these problems can be eliminated with moderate extra computational effort. Here 
we present the pragmatic solution of using larger R c such that the remaining error due 
to the presence of the Coulomb potential becomes acceptable. 

Fig. 3 shows a spectrum calculated for the Hydrogen atom with a FWHM T = 
20 optical cycle pulse at 800 nm wave length and peak intensity 10 14 W / / cm 2 . All 
discretization errors can be controlled in the same way as discussed for the short range 
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Figure 3. Photo-electron energy spectra for the hydrogen atom (upper panel) 
obtained with surface radius R c = 110 (black) and R c = 140 (red). Middle panel: 
error estimate by comparing the spectra according to Eq. (22). Lower panels: blow-up 
of the spectra from the upper panel on a linear scale. Pulse parameters: A = 800 nm, 
T = 20 optical cycles, intensity = lQ^W/cm 2 . 
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potential. The error is dominated by the dependence on the surface radius R c : on an 
absolute (logarithmic) scale, two calculations with R c — 110 and R c = 140 are hardly 
discernable. The error level of the calculation with R c = 110 and 180 discretization 
points is < 10% and decreases slowly as R c increases. In the linear plot of the spectra 
(lowest panels of Fig. 3) we see that the largest errors occur at lower energies due to the 
larger influence of the weak Coulomb tail on low energy scattering states. Agreement 
at intermediate energies is near perfect. The increase of relative errors at the highest 
energies is due to a displacement of peaks caused slightly incorrect dispersion due to 
spatial discretization. 

Fig. 4 shows angle-resolved photo-electon sprectra for FWHM durations of T=10, 
20 and 30 optical cycles. In the region up to 10 U p ~ 60 eV circular structures with their 
centers offset along the polarization axis are clearly distinguishable, best visible at the 
shortes pulse T = 10. These structures were first explained in Ref. [17] and are due to 
re-scattering. The intensity around each circle is related to the electron-ion scattering 
differential cross section [18]. In the enlarged plots of the energy region up to energies 
of 2 U p 12 eV we reproduce rings and fan-like structures that are related to particular 
partial waves, as has been extensively discussed in Ref. [19]. Sub-structures between 
the photo-electron energy peaks are clearly resolved. These are related to the pulse 
envelope: the spacing decreases and number of peaks increase as the pulse-duration 
increases. In different wording one can say that they are caused by interference between 
rising and trailing edges of the pulse [20] . 

For comparison with results reported in [9], we also computed the spectrum at 
the lower intensity of 5 x 10 13 W/cm 2 and T = 10 cycles FWHM (corresponding 20 
cycles total pulse duration). Fig. 5 shows our result for the total photo-ionization 
spectrum obtained with only 192 discretization points and 25 angular momenta. We 
have estimated the accuracy of our results to be < 10% throughout the spectrum 
using the same procedures as for Fig. 3. Our result qualitatively differs from Fig. 2 
in Ref. [9], where a surprsing irregularity appears near 10 U p , while no such structure 
exists in our sprectrum. Unfortunately no detailed discussion of accuracy or convergence 
is given in [9]. One possible source of the discrepancy is insufficient discretization. With 
peak energy of 40 U p ~ 4a. u. included in the calculation the maximal momenta are 
Pmax ~ 2.8 cl.u.. During 20 optical cycles (2200 a.u. of time), these electrons move 
to distances ~ 6000 a. u.. With the 2000 discretization points used in [9] one gets 
average momentum grid spacing of Ap = p max /2000 and an effective spatial box size 
~ 2n/Ap PS 4500, which appears somewhat below the necessary limit. However, the 
uneven distribution of grid points and additional spectral cuts in the energy domain used 
in [9] make it difficult to carry this analysis further. Rather, a systematic convergence 
study would be required. The present number of discretization points also compares 
well to the 1000 ~ 2000 discretization coefficients used in [11] at photon energies 
0.3 ~ 0.7 a.u.. The results of the benchmark calculations Refs. [21] obtained with 
box size of 3000 atomic units and about 4000 discretization points at somewhat shorter 
wave length of 620 nm could be reproduced using only 200 radial discretization point 




Figure 4. Angle resolved photo-electron energy spectra for the hydrogen atom at 
A = 800 nm, intensity 10 U W/ cm 2 and FWHM pulse durations T = 10, 20, and 30 
(upper panels). In the lower panels, the region up to 2U P ~ 12 eV is enlarged. 



t-SURFF: photo- electron spectra from minimal volumes 



13 




1 







1 







1 







0.4 0.8 1.2 

Photoelectron energy (a.u.) 



1.6 



Figure 5. Photo-electron energy spectrum for a FWHM T = 10 pulse at 800 nm and 
peak intensity I = 5 X 10 13 W/cm 2 . Results obtained with 192 radial discretization 
points are accurate to a few percent and never exceed 10% in the whole range shown. 
The pulse parameters agree with those used for Fig. 2 in Ref. [9] . Near 10 U p there 
apears a striking qualitative difference between our result and Ref. [9]. (See text for a 
discussion.) 



up to energies of 10 U p , except at the lowest energies, where our method is limited by 
use of Volkov solutions beyond the surface radius R c < 150a. u.. 

6. Extensions of the method 

6.1. Handling long-range potentials 

In the section above we were able to obtain good results for quite demanding laser 
parameters using only R c ph 140. Still, the problem remains big and in particular 
when considering extension of the approach to multi-electron systems, mult i- dimensional 
simulation volumes would quickly exhaust computational resources. In turn, a reduction 
of R c to the necessary minimum of 20 ~ 50 set by the quiver radius, would constitute an 
essential gain, possibly deciding about the feasibility of the multi-electron calculation. 
Here we show how to correctly handle long-range potentials in t-SURFF. 

For computing exact spectra, we must know the exact solution of the TDSE beyond 
R c . In the derivation of Eq. (15) we have only used H(t) = H c (t) for \r\ > R c and that 
we can by some means obtain accurate solutions Xk(r,t) for the TDSE with H c (t). A 
suitable H c is obtained by using in Eq. (19) the potential 




1/rforr > R c . 



(23) 



t-SURFF: photo -electron spectra from minimal volumes 14 

Partial wave solutions (j>k,i(f) for the field-free case A = with V c are the spherical Bessel 
functions up to R c which are connected smoothly to values and derivatives of regular 
and irregular Coulomb functions in the region r > R c , where one must pay attention 
to proper (^-normalization. With non-zero field S(t), no exact solution is knwon. We 
found that a simple Coulomb- Volkov approximation [22] for the time-dependence is 
insufficient: we could not observe any acceleration of convergence by replacing the plane 
waves of the Volkov solution with the field-free scattering solutions. Lacking a reliable 
approximation for the scattering solution, we must solve the TDSE for the x%{r, t) with 
the final condition 

XkAr,T) = Mr)- (24) 

As the potential is weak, little actual scattering occurs and accurate solutions can be 
obtained by expanding into <f>k',i'(r) for small intervals around the asymptotic radial and 
angular momenta, k' G [k — Ak, k + Ak] and I' e [I — A/, I + A/]. Numerical results 
using this procedure will be presented elsewhere. 

Finally there remains the poblem that Rydberg states may become occupied, which 
have non-neglible amplitude at r = R c . The effect of bound states on the integrals (15) 
is slow, oscillatory convergence to the asymptotic value at T — > oo. Again, there is 
an efficient and rather pragmatic solution to this problem: by averaging the value over 
several optical cycles, we find that convergence is speeded up and the reported accuracies 
are reached quickly. 

If the Rydberg states are known exactly, we can remove them after the pulse is 
over. This removes all oscillations and the asymptotic value is reached rapidly. For 
the procedure we define the projector onto the (field free) Rydberg states \n) and its 
complement as 

P := ^2 \n)(n\ and Q := 1- P. (25) 

n 

A simple calculation shows that the spectral amplitude with the Rydberg states removed 
is 

~ " ') + 



(xstT)\o(Rc)Q\*(T)) = i [ ° dt(xM[ H MO(R c )]Mt))- 

J ~ oo 

[ T dt( Xk -\[H c (t),e(R c )}QMt)) - ( Xj ;(T )\e(R c )PMT )), (26) 



where T is any time after the end of the pulse. If high precision is wanted and if long 
time-propagation is costly, the extra effort of implementing the explicit projection (26) 
may be justified. 

6.2. Single-ionization of multi- electron systems 



The procedure can be extended to describe single-ionization of multi-electron systems. 
The main new feature here is that we have several ionization channels, depending on 
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the state in which the ion is left behind. The spectral density in a specific channel c can 
be written, as before, as 

^(k) = \(x crk (T)\9(R c )\UT))\ 2 . (27) 
The asymptotic channel wave function x c fullfills the channel TDSE 

*T t Xcjfr) = H e (t) X Jt), (28) 

with the channel Hamiltonian 

H c (t) = ^[-iV - A(t)} 2 ® H ion (t) (29) 

where for simplicity we neglect the ionic Coulomb potential. Note that we do not need 
to explicitly anti-symmetrize Xc if is anti-symmetric. A general solution of (28) has 
the form 

XcS {t) = (27r)- 3 /V^V fe> ® c (t), (30) 

where </> c (t) solves the ionic TDSE with Hamiltonian Hi on (t) and a final state condition 
that specifies the ionic state c of the channel: 

H lon (T)<p c (T) = E c <t> c {T). (31) 

Assuming that double-ionization is negligible, the further steps for computing the 
spectral amplitude as a time-integral are the same as in the single-electron case. In 
addition to ^ s (t) we must also compute c (t). Values and derivatives of the c-channel 
scattering wave function tp c (r,t) := (<f) c (t)\^ s (t)) at \r\ = R c can be stored for a 
sufficiently dense time-grid. Because of the stronger binding of electrons in the ion, 
computing the ionic wave function <f> c (t) usually requires much xsless effort than obtaing 
*.(*)■ 



6.3. Double-ionization 

For double-ionzation spectra, we must know the two-electron solution in the asymptotic 
region. This may be less difficult than what it appears at first glance. Neglecting 
the ionic potentials, we have a Volkov solution for the center of mass coordinate and 
Coulomb waves for the relative coordinate. What is left to do is to match that solution 
with an accurate solution on a five-dimensional surface, where all reflections from the 
simulation box boundaries are carefully suppressed. While this procedure has not been 
worked out in detail, it may be well feasible and bring IR two-electron spectra from 
realm of herioc super-large scale calculations [5] to managable size, allowing systematic 
studies. 



7. Conclusions and outlook 



We have shown that the unfavorable scaling for the computation of photo-electron 
spectra with laser wave length can be largely overcome with t-SURFF, which picks 
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up the exact solution at some finite surface and beyond that surface exploits knowledge 
about the long-range behavior of the solutions of the TDSE. The reduction of problem 
size is particularly striking for atomic binding potentials with finite range, when the 
Volkov-solutions become exact at distances where the potential is zero. In that case, 
the box sizes can be reduced to approximately the range of the potential plus the electron 
quiver amplitude. Electrons that move beyond that range will never scatter and will 
exactly follow the Volkov solution. For our parameters, the wave function expands to 
several thousand atomic units during the pulse and correspondingly large boxes would 
be needed, if a spectral analysis of the wave function were performed after the end of 
the pulse. In contrast, we could present < 1% accurate spectra up to energies of 120 eV 
using a box size of only about 30 atomic units and as little as 75 discretization points 
per partial wave. With only a few more points, much higher accurcacies can be reached. 

Instrumental for the application is the traceless absorption of the wave function 
beyond the surface, which is provided by the irECS method introduced in a preceding 
publication [16]. The good performance of irECS for one-dimensional wave functions and 
for high-harmonic signals from three-dimensional calculations presented in [16] could be 
confirmed also for the much more delicate observable of angle-resolved photo-electron 
spectra. 

For the long-range Coulomb potential, the Volkov solutions are not exact 
asymptotically, let alone at any finite distance. In physical language, the electron will 
scatter in the long-range tail of the Coulomb potential and pick up more energy from the 
laser field and therefore we cannot predict its final energy before the pulse is over. An 
attempt to approximate the asymptotic behavior by Coulomb- Volkov instead of the pure 
Volkov solutions was futile. In a pragmatic approach, we could show that with a surface 
at distances of R c = 100 ~ 140 atomic units, still impressively accurate spectra can be 
obtained using Volkov solutions as approximation to the exact asymptotic solutions. 

For single-electron systems and with moderate accuracy requirements, simulation 
volumes on the scale of ~ 100 atomic units are quite acceptable. For experimentally 
interesting multi-electron systems, we have proposed to further reduce box-sizes by 
solving the asymptotic laser-assisted scattering problem numerically. This may be done 
efficiently for an asymptotic Hamiltonian that only includes scattering at distances > R c 
and dismisses the main part of scattering from near the Coulomb singularity. The need 
to solve this weak scattering problem enhances the complexity of coding and significantly 
increases computation times. However, for few-electron systems, this extra effort is far 
compensated by an expected reduction to box sizes to as little as 20 ~ 50 atomic units. 
We have also formulated the extension of t-SURFF to single- and double-ionization of 
multi-electron systems. A numerical demonstration of these methods will be the subject 
of future work. 

In summary, t-SURFF, while producing highly accurate results, can reduce the box 
size for computing IR photo-electron spectra for single electron systems by one order 
of magnitude or more. This drastic reduction of box-sizes is particularly important for 
very long laser wave length and for multi-electron systems. For systems with two and 



t-SURFF: photo- electron spectra from minimal volumes 



17 



more electrons, we believe, it opens a root to computing accurate IR photo-electron 
momentum spectra. 
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